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Abstract 

We investigate the effects of the nuclear equation of state (EoS) to the neutron star cooling. New 
era for nuclear EoS has begun after the discovery of ~ 2Mq neutron stars PSR J1614—2230 
and PSR J0348-I-0432 [U |2]. Also recent works on the mass and radius of neutron stars from 
low-mass X-ray binaries [3] strongly constrain the EoS of nuclear matter. On the other hand, 
observations of the neutron star in Cassiopeia A (Cas A) more than 10 years confirmed the existence 
of nuclear superfluidity Nuclear superfluidity reduces the heat capacities as well as neutrino 

emissivities. With nuclear superfluidity the neutrino emission processes are highly suppressed, and 
the existence of superfluidity makes the cooling path quite different from that of the standard 
cooling process. Superfluidity also allows new neutrino emission process, which is called ‘Pair 
Breaking and Formation’ (PBF). PBF is a fast cooling process and can explain the fast cooling rate 
of neutron star in Cas A. Therefore, it is essential to add the superfluidity effect in the neutron star 
cooling process. In this work, we simulate neutron star cooling curves using both non-relativistic 
and relativistic nuclear models. The existence of too early direct Urea process shows that some of 
nuclear models do not fit for the cooling simulation. After this first selection process, the nuclear 
pairing gaps are searched using the observational neutron star’s age and temperature data. 
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1. Introduction 

A neutron star (NS) is born as a result of the core collapse supernova explosion if the initial 
mass of the main sequence star is around 8 to 20 solar mass (Mq). Various EoS, relations between 
the density and the pressure, are used to describe the interior of neutron stars. The resulting central 
density of a neutron star is expected to reach up to several times of the nuclear saturation density 
(no — 0.16 fm“^). Hence, a neutron star is one of the best astrophysical laboratories to study the 
physics of the extremely dense nuclear matter. 

Figure shows the structure of neutron stars in the theoretical point of view. Outside the 
envelope, there exists a very thin atmosphere which is composed of hydrogen and heavy elements 
[6]. Depending on the temperature, its thickness ranges from a few centimeters to a few millimeters 
and the thermal radiation is expected to occur in this region. The envelope plays a role as a thermal 
insulator [7j bewteen the outer crust and the atmosphere. In the outer crust, heavy nuclei exist 
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Figure 1: (Color online) A schematic of neutron star structure. 


with free gas of electrons, and the BCC (Body Centered Cubic) structure of nuclei is expected. 
As the density increases, neutrons drip out of heavy nuclei and can exist as free gas in the inner 
crust. The neutron superfluidity can also exist in the inner crust. In the boundary of the 
outer core and the inner crust, nuclear pasta phase may appear due to the nature of Coulomb 
interaction 019 ], which favours a special geometrical shape to reduce the free energy. Above half 
of the nuclear saturation density uniform nuclear matter exists. Outer and inner cores occupy 
roughly 70 % of neutron star’s volume. In the outer core the uniform nuclear matter can exist and 
neutron ^P 2 and proton superfluidities are expected. In the inner core exotic states with new 
degrees of freedom, such as kaons, hyperons, and deconfined quarks might appear. However, the 
inner structure of neutron star is still unknown due to the theoretical uncertainties, especially at 
the densities far beyond the nucleon saturation density. 

The mass distributions of observed neutron stars are summarized in Table 1 of Lattimer m- 
Note that the mass distribution of NS-NS binaries has very sharp peaks since the error bars of mass 
measurements are relatively small, and the well measured neutron stars in NS-NS binaries have 
masses < 1.5Mq. On the other hand, 2Mq neutron stars have been observed in NS-WD (white 
dwarf) binaries. This implies that the mass distribution may depend on the binary evolution in 
addition to the neutron star equation of states m- Recent observations of ~ 2Mq neutron stars 
ruled out many soft equations of states with which the maximum mass of neutron star becomes less 
than 2Mq mm- Another constraint on the nuclear equation of state comes from the X-ray binary 
observations and their analysis 0 . Mass distribution in X-ray binaries is very broad without sharp 
peak indicating that the uncertainties in the mass estimations are very large. 

In Table and Fig. we plot the surface temperatures (T“) and photon luminosities {L°°) 
of 18 observed neutron stars. Note that three different models are used to estimate the effective 
surface temperatures of neutron stars at infinity. For black body (BB) model, the photon luminosity 
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{L°°) and the effective surface temperature (T^) of a neutron star at infinity are related by 

L“ = 47r(i?“)VsB(rr)" (1) 

where is the radius seen by an observer at infinity and ctsb is the Stefan-Boltzmann constant. 
We can determine R°° by inserting and L°° in the Table into Eq. 0- 

In recent works na 13], we have investigated the behavior of nuclear EoS at high densities 
by calcnlating the mass and radius of neutron stars with several Skyrme force models. We could 
confirm that the models that satisfy the large mass observation mm are consistent with the mass- 
radius zone in [3]. Moreover, the conclusion doesn’t change even if we include exotic degrees of 
freedom such as kaon condensation |12] or hyperons m- On the other hand, nuclear EoS is one of 
the key ingredients that determine the thermal evolution of neutron stars. Available nuclear models 
predict very diverse mass-radius relations, so it is well expected that the cooling behavior will be 
sensitive to nuclear models. Therefore, cooling of nentron stars provides a multi-test ground for 
the nuclear models; whether a model good for mass-radius relation is consistent with the observed 
temperature data or not. 

In this work, we simulate the cooling curve of neutron stars with varions EoS obtained from 
both relativistic and non-relativistic models. Without the superfluidity effect, the direct Urea 
process is a good indicator whether a specific nuclear model is suitable for neutron star cooling 
simulation. If the temperature drop is too fast even for the low mass neutron stars ( 1.2 Mq), 

such EoS is not consistent with most of the neutron star temperature observations. Even though 
the existence of superfluidity highly suppresses neutrino emissivity, once the direct Urea process 
is turned on before the critical temperature, the superfluidity cannot slow down the temperature 
drop so that the temperature of neutron stars is much below any observational data.Thus, the 
study of cooling curve of neutron star can provide hints about the inner structure of neutron stars 
and the EoS of dense nuclear matter. In addition to the nuclear EoS, physical conditions snch as 
the composition of elements in the envelope, and existence of superfluidity in the core play crucial 
roles in determining the cooling curve of neutorn stars. In general, chemical composition of the 
envelope seldom affects the cooling mechanism, but the surface temperature depends sensitively on 
whether the elements in the envelope are light or heavy. On the other hand, superfluidity directly 
determines the cooling rate. If nucleons form a cooper pair and transit to a snperfluidic state, the 
rate of neutrino emissivity is suppressed exponentially. This may lead to a very slow cooling rate. 
However, below a critical temperature, creation and destruction of cooper pairs ignite a fast cooling 
mechanism, and this can give an abrupt decrease of temperature. Recent literature succeeds to 
reproduce the cooling curve of Gas A in terms of PBF HE]. In this work, we incorporate PBF to 
various nuclear models and explore the extent to which the models can reconcile with PBF . 

The structure of the paper is as follow. In Section we summarize the equation of state of 
nuclear matter which we used in the study of thermal evolution of neutron stars. In this work, 
we do not consider exotic matter or quark matter in the core of neutron stars but assume that 
the core is composed of nucleons (protons and neutrons) and leptons (electrons and muons) in 
the form of uniform matter. In Section we present the basic equations for neutron star cooling. 
Simple critical temperature formula is proposed and we analyze the pairing effects to neutron star 
cooling. In Section]^ we present the cooling curves for the standard cooling in which superfluidity 
is neglected with various nuclear models. In Section we discuss the effect of the superfluidity to 
the cooling process. In Section we give conclusions from neutron star cooling enrves combined 
with varios EoS. In Appendix, we summarize numerical methods to solve the diffusion equations 
for the thermal evolution of neutron stars. 
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Table 1: Thermal emission from isolated neutron stars. Temperatures were obtained using three 
different models; hydrogen atmosphere (HA), magnetized hydrogen atmosphere (mHA), and black- 
body (BB) model. For sources with no. 2 ~ 7, two different models were used to link the effective 
temperature and the photon luminosity. Sources with no. 16 ~ 17 have the limited observational 
data, thus have only one limit. 


No 

Source 

Log(Gd/yr) 

Log(4m/yr) 

Log(Tr/K) 

Log(L°°/erg s °) 

Model 

Ref. 

1 

PSR 11119-6127 

3.20 

- 

o 

00 

1 + 
o o 
b b 

32.88 - 33.66 

mHA 

m 

2 

RX 10822-42471' 

3.90 

o r7+0.04 

0.0 i _o .04 

(4 24+®'^'^ 
O.Z^_0.04 

a f;c:+o.04 
^■^^-0.04 

33.85 - 34.00 

33.60 - 33.90 

HA 

BB 

[mile] 

3 

IE 1207.4-5209 

r ro + 0.44 
^•^'^-0.19 

Q qix+ 0.48 
O.00_o.48 

6 - 2 ilS);!)^ 

6.481°;°) 

33.27 - 33.74 

33.27 - 33.74 

HA 

BB 

[niiis] 

4 

PSR 11357-6429 

3.86 

- 

5-88iro4 

a. oq+0.05 

U.ZO_Q Q5 

32.46 - 32.80 

32.35 - 32.76 

mHA 

BB 

m 

5 

RX 10002-H6246 

- 

O.yu_o.o 8 

a nQ+0.03 
'^■'^'=^-0.03 

« 1 c:+ 0.11 
O. 10 _o.ii 

33.08 - 33.33 

32.18 - 32.81 

HA 

BB 

[mile] 

6 

PSR B0833-451' 

4.05 


K qq+ 0.02 
^'^'^- 0.02 

6.18l°;°02 

32.41 - 32.70 

32.04 - 32.32 

mHA 

BB 

nnnn] 

7 

PSR B1706-44 

4.24 

- 

5.801°;)° 

no + 0.04 
'^■"^^-0.04 

31.81 - 32.93 

32.48 - 33.08 

mHA 

BB 

mnn] 

8 

PSR 10538-t2817 

4 47+0-05 
+•+'-0.06 

- 

r q^+0.08 

o.y^_o.o 8 

32.32 - 33.33 

mHA 

m 

9 

PSR B2334-t61 

4.61 

- 

^■^+-0.08 

31.93 - 32.96 

mHA 

m 

10 

PSR B0656-tl4 

5.04 

- 

e; 71 +0.03 
^ -*--0.04 

32.18 - 32.97 

BB 

[241 [I 6 ] 

11 

PSR B0633-tl748l' 

5.53 

- 

e; 7 C+O.O 4 
^^-0.05 

30.85 - 31.51 

BB 

113 dS] 

12 

RX 11856.4-3754 

- 

r yQ + 0.05 
'^-0.25 

5-63irol 

31.32 - 32.35 

mHA 

nans 

13 

PSR B1055-52 

5.73 

- 

5-88iro8 

32.05 - 33.08 

BB 

HZldS] 

14 

PSR 10243-t2740 

6.08 

- 

r ^^+0.08 
^■^+-0.08 

29.10 - 30.13 

mHA 

m 

15 

RX 10720.4-3125 

6.11 

- 

■t 70+0-08 

0 . m_o.o 8 

31.37 - 32.40 

HA 

m 

16 

PSR 10205-t6449H 

- 

2.91 

< 6.01 

< 33.29 

BB 

m 

17 

PSR B0531-t2lV 

- 

3.0 

< 6.30 

< 34.45 

BB 

m 

18 

RX 10007.0-^730311 

- 

4.0-4.4 

< 5.82 

< 32.54 

BB 

m 


t Alternative names: Puppis A (PSR B0822-4247), Vela (PSR B0833-45), Geminga (PSR B0633-H1748). 

PSR J0205-I-6449 is a pulsar in supernova remnant 3C 58, PSR B0531-I-21 is in SN 1054 in Crab Nebula, 
and RX 10007.0-^7303 is in the CTAl. 
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Figure 2: (Color online) The effective temperature and the photon luminosity for observed neutron 
stars. Left panel : the effective temperature at infinity. Right panel : photon luminosity seen by 
the observer at infinity. Note that three different models are used to link the effective temperature 
and the luminosity as in Table 


2. Neutron Star Equation of State 

Nuclear matter properties beyond nuclear saturation density are still unknown and different 
EoSs give quite different masses and radii of neutron stars. In order to investigate the properties of 
neutron star matter, we first consider both relativistic and non-relativistic models for the neutron 
star core, which are consistent with 2.0Mq neutron stars mm- We consider the crust of neutron 
star separately because heavy nuclei can exist in the crust. The properties of neutron star crust is 
very important in understanding the neutron star properties in low-mass X-ray binaries. 


2.1. Non-relativistic nuclear force model : energy density functional 

For the non-relativistic nuclear force model, we use Skyrme force model to obtain the equation 
of state for nuclear matter [32]. In Skyrme force model, the interaction between two nucleons has 
the form of 


vsF{ri,rj) = to{l -\- xoPa)S{ri - Tj) + ^-^{1 + xiPa) ()(ri - rj)k^k'f^(5(rj - Fj) 

+t2{l + X2Pa)k'^ • (5(ri - rj)k -k + X3pa)n°‘S{ri - rj) 

+iWQiP5{ri - Tj) X k • (d-j -k a-j ), 


( 2 ) 


where P^ = |(1 -k • <Tj) is the spin-exchange operator, tj, Xi and a are the parameters of the 
interactions, and k is defined as 

k=^(V*-V,). (3) 

Note that the interaction contains terms up to quadratic in derivatives, ta-term is added to include 
many body effect beyond quadratic in density n, and kFo-term gives the spin-orbit interaction 
which is important to explain the nuclear structure. Parameters used for Skyrme force models are 
summarized in Table [2j 
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Table 2; Parameters used for Skyrme force models. 



SLy4 

SkI4 

SGI 

SV 

TOV min 

LS 220 

to (MeV- fm^) 

-2488.913 

-1855.83 

-1603.0 

-1248.3 

-2129.735 

-2182.73 

ti (MeV- fm®) 

486.818 

473.829 

515.9 

970.6 

305.398 

0 

t 2 (MeV- fm^) 

-546.395 

1006.855 

84.5 

107.22 

362.532 

0 

ts (MeV- fm3+3“) 

13777.0 

9703.607 

8000.0 

0.0 

13957.064 

14960.296 

Xq 

0.8340 

0.405082 

-0.02 

-0.17 

0.169949 

-0.3121 

Xl 

-0.3438 

-2.889148 

-0.5 

0.0 

-3.399480 

0 

X2 

- 1.0 

-1.325150 

-1.731 

0.0 

-1.782177 

0 

X3 

1.3540 

1.145203 

0.1381 

1.0 

0.402634 

-0.5 

a 

1/6 

0.25 

1/3 

1 

0.2504 

0.2602 

h? 

2mr) 

20.7525 

20.7525 

20.7525 

20.7525 

20.7498207 

20.7484 

V 

2m„ 

20.7525 

20.7525 

20.7525 

20.7525 

20.7212601 

20.7484 


In Hartree-Fock level, the total energy can be expressed as 

E = ^{i\i\j)Pji + ^ X] ^ijkiPkiPij (4) 

ij ijkl 

where t is the kinetic energy operator and 

Vijkl = - PaPrPT)\kl). (5) 

Here Pr is the parity operator and Pr = ^{1 + n ■ tj) is the iso-spin exchange operator. Due to the 
zero range property of Skyrme force model, the total energy can be easily obtained as 

E = J d\£ = I Pr{£B+£c + £9+£j) ( 6 ) 

where Ts is the bulk part contribution, £c is the Coulomb contribution, £g is the contribution 
from the density gradient term, and £j is the contribution from the spin-orbit term. For a uniform 
matter in the neutron star core, Tb is dominant. Hence the energy density can be approximated 

as [Mj 
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Table 3: Nuclear matter properties at the saturation density (no). Upper 6 models correspond 
to non-relativistic Skyrme force models and lower 6 models correspond to relativistic mean field 
models. B is the binding energy of the symmetric nuclear matter, is the symmetry energy, L 
is the slope of the symmetry energy, K is the compression modulus, and is Landau effective 
nucleon mass (effective chemical potential). 


Model 

no (fm ^) 

B (MeV) 

Sy (MeV) 

L (MeV) 

K (MeV) 


Ref 

SLy4 

0.160 

16.0 

32.0 

45.9 

230 

0.694 

m 

SkI4 

0.160 

16.0 

29.5 

60.4 

248 

0.649 

136] 

SGI 

0.155 

15.9 

28.3 

63.9 

262 

0.608 

137| 

sv 

0.155 

16.1 

32.8 

96.1 

306 

0.383 

m 

TOV-min 

0.161 

15.9 

32.3 

76.2 

222 

0.934 

[39| 

LS220 

0.155 

16.0 

28.6 

73.1 

220 

1.000 

m 

lU FSU 

0.155 

16.4 

31.3 

47.2 

231 

0.669 

m 

DD-MEd 

0.152 

16.1 

32.4 

52.9 

219 

0.668 

m 

SFHo 

0.158 

16.2 

31.6 

47.1 

245 

0.810 

m 

NLp 

0.160 

16.1 

30.4 

84.6 

241 

0.800 

m 

TMA 

0.147 

16.0 

30.7 

90.1 

318 

0.691 

m 

NL3 

0.148 

16.2 

37.3 

118 

272 

0.655 

m 


where and nip are neutron and proton masses, Un and Up are the densities of neutrons and 
protons, the total baryon density n = + np, and and Tp are kinetic energy densities of 

neutrons and protons, respectively. The pressure can be obtained by taking a density derivative of 
energy per baryon, 


P = n 


• d{£/n) 

dn 


( 8 ) 


In the upper part of Table we summarize the basic properties of nuclear matter for Skyrme force 
models which are used in this work. In the upper part of Figure we plot the pressure for Skyrme 
force models for both symmetric and pure neutron matter. In the left panel of Figure masses 
and radii of neutron stars are summarized for Skyrme force models. 


2.2. Relativistic mean field model 

Relativistic mean field models have been very successful in explaining the nuclear properties, 
such as binding energy, density profile, root mean square radius, etc. Based on these successes, 
many efforts have been given to build RMF models which are suitable for both finite nuclei and 
bulk nuclear matter properties. The typical RMF models are described by the Lagrangian |34j . 
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(9) 
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Figure 3: (Color online) Pressure of symmetric nuclear matter and pure neutron matter from 
both non-relativistic Skyrme force models (upper panels) and relativistic mean field models (lower 
panels). The shaded area is the result from Ref. |47) . For the pure neutron matter (right panels), 
the upper (lower) shaded area in each plot represents the stiff (soft) equation of state. 





























Figure 4: (Color online) Neutron star’s mass and radius relation from non-relativistic Skyrme force 
models (left) and relativistic mean field models (right). Thick horizontal lines indicate the masses 
of PSR J1614—2230 and PSR J0348+0432 [U [2]. The shaded area is the most probable mass and 
radius, Icr and 2a region, from the analysis of Steiner et al.^. 


where a is the scalar field, is the vector-isoscalar field, is the vector-isovector field, A is 
the photon field, 6 is the scalar-isovector field, — di,A^, = d^ojy — Rfiu = 

d^by — dyb^^ and 14ff is the general effective potential for meson fields. N The equations of motion 
for meson fields can be obtained using the Euler-Lagrange equation 


3 ( \ - 

^\d{d^(l))) d(f) 

where cj) = a,oj^,bf,6i. By taking expectation value of each field, one can define the scaled meson 
mean fields <1> = ga{c’'), W = gi^{uj^), and R = gp{b^)- Then the equations of motion for the uniform 
nuclear matter become 



Us 


n 



1 dVes{^,W,R) 

-I 54> 

1 dVes{^,W,R) 

^2 dw 

1 ^ dV,s{^,W,R) 


( 11 ) 

( 12 ) 

(13) 


where the scaled coupling Cj = gijmi, the baryon scalar density Us = {‘4’^), baryon density n = 
-)- kp^)/{3TT‘^), and baryon isovector density ns = = {kp^ — kp^)/{37r‘^). The 

pressure and energy density for nuclear matter can be obtained as 


P 



J_$2 , + —R2 

2c2 +202"^ +2c2" 


Ves{^,W,R), 


^In some literature, meson mass terms are also included in the effective potential. However, in this work, mass 
terms are explicitly specified and the effective potentials have only higher order interaction terms beyond mass term. 
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1 


+Wn+-Rp3, 


(14) 


V;fr('7,w'"w^,6'"-6^) = ^{gaaf+:^{g^a)'^-^g^{uj^uj^f-^gp{P-bf,f-f{a,u)>^uii,)glb^^-b^, (15) 


(16) 


where m*j^ = mN — gaO’^gs^s (+ : proton, — : neutron). Various forms of meson effective potentials 
14ft(4>, W, R) are used in the literature, 

with 

6 3 

/(a, w^w^) = OiO-* + bj{(jj^ujij,y . 

i=l j=l 

The values of the parameters in various models are summarized in Table In the density dependent 
coupling constant model (DD-ME(5), the coupling constant has the form of 

9\{n) = g\{no)sx{x ), (17) 

where x = n/no, X = a, ui, p, 5 and 

1 + 6a (x + dx) 


sx{x) = ax: 


(18) 


T + CAi(x + eA) ' 

Numerical values of parameters are summarized in Table Note that gx^s in Table are equal to 
gx{no) in Eq. ( ]17[ ). In the lower part of Tablewe summarize the basic nuclear matter properties 
of RMF models selected in this work. In the lower part of Figure we plot the pressure for various 
RMF models for both symmetric and pure neutron matter. In the right panel of Figure]^ masses 
and radii of neutron stars are summarized for RMF models. 


2.3. Neutron star crust 

In the crust of neutron stars, heavy nuclei are expected to exist with free gas of neutrons and 
electrons. A simple but appropriate description of this situation is completed using liquid droplet 
formalism nano]. The total energy density (without electron contribution) is given by 


F = UUifi + + PsT^n] + ^(rArniXie)^c(u) + (1 - u)nnofo , 

rN 5 


(19) 


where u is the volume fraction of heavy nuclei to Wigner-Seitz cell, n* is the density inside of heavy 
nuclei, /* is the energy per baryon of the heavy nuclei, s{u) is the surface shape factor, r^f is the 
radius of heavy nuclei, a{x) is a surface tension as a function of proton fraction x, ps is the neutron 
chemical potential on the surface, Un is the areal neutron density on the surface, x, is the proton 
fraction of heavy nuclei, c{u) is the Coulomb shape function, n„o is the neutron density outside 
of heavy nuclei, and fo is the energy per baryon outside of the heavy nuclei. Minimizing energy 
density, we have four equations to solve 


uuiXi — nYp = 0 , 

2/3 

uni + X—Ti't'n + (1 - u)nno - n = 0, 

0(7 

pni Pno — 0, 


(20) 
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Table 4: Paramters for RMF models. We multiply gp by the factor of 2 in DD-ME(5 ,NL/9, and 
TMA models since the definition of Lagrangian is different from the standard Lagrangian here we 
used. For DD-ME5, all the quoted coupling constants {ga, Qu, Qp, ge) are the ones calculated at the 
saturation density. 


Par am. 

Unit 

lU-FSU 

DD-MES 

SHFo 

NLp 

TMA 

NL3 

rria 

MeV 

491.500 

556.158 

457.286 

508.194 

519.151 

508.194 

rriuj 

MeV 

782.500 

783.000 

762.500 

782.501 

781.950 

782.501 

nip 

MeV 

763.000 

763.000 

770.000 

763.000 

768.100 

763.000 

ms 

MeV 

- 

983.000 

- 

- 

- 

- 

9<j 


9.9713 

10.3325 

7.53631 

8.27739 

10.055 

10.217 

9iJj 


13.0321 

12.2904 

8.78166 

9.23205 

12.842 

12.868 

9p 


13.5900 

2x6.3128 

9.38351 

2x3.76877 

2x3.800 

8.92188 

95 


- 

7.1520 

- 

- 

- 

- 

K 

fm“^ 

0.017133 

- 

0.0710487 

0.066 

6.4529 X lO""^ 

3.8599 

A 


0.000296 

- 

0.0245322 

-0.0288 

0.0228112 

-0.015905 

c 


0.03 

- 

-0.0017013 

- 

0.0334419 

- 

e 


- 

- 

0.0034525 

- 

- 

- 

ai 

fm“^ 

- 

- 

-0.23016 

- 

- 

- 

02 


- 

- 

0.57972 

- 

- 

- 

03 

fm 

- 

- 

0.34446 

- 

- 

- 

04 

fm^ 

- 

- 

3.4593 

- 

- 

- 

05 

fm^ 

- 

- 

1.3473 

- 

- 

- 

06 

fm^ 

- 

- 

0.66061 

- 

- 

- 

hi 


7.81241 

- 

5.8729 

- 

- 

- 

h2 

fm^ 

- 

- 

-1.6442 

- 

- 

- 

bs 

fm^ 

- 

- 

314.64 

- 

- 

- 


Table 5: DD-ME(i density dependent parameters [12]. A indicates type of meson. 


A 

ax 

bx 

cx 

dx 

ex 

a 

1.3927 

0.1901 

0.3679 

0.9519 

0.9519 

u 

1.4089 

0.1698 

0.3429 

0.9860 

0.9860 

P 

1.8877 

0.0651 

0.3469 

0.9417 

0.9737 

5 

1.5178 

0.3262 

0.6041 

0.4257 

0.5885 
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with four unknowns, u, n*, Uno, and Xj. Here, /3 = 9(7re^x?n?(T^/15)^/^, T>' = dVjdu, and V = 
[c(tt)s^(u)]^/^ is a geometric shape function which corresponds to nuclear pasta phase in liquid 
droplet model [8lll0]. Pi (Pq) is pressure inside (outside) of the heavy nuclei and the total pressure 
is given by 

P = P^-p{p- uV'). ( 21 ) 

The boundary between the crust and the core can be found by comparing the energy density or 
energy per baryon of uniform nuclear matter and the heavy nuclei with free neutron and proton 
gas. In general case, the energy difference between two phases near the boundary is so small that 
the pressure difference is negligible [9j. 


3. Neutron Star Cooling Mechanism 


Thermal evolution of neutron star can be obtained by solving the diffusion equations 


L'f' 

1 




rc^ 


dr 




rc^ 


dr 


q dT 

e'^3 dt ’ 


( 22 ) 

(23) 


where is the local photon luminosity, T is the local temperature, m = m{r) is the enclosed mass, 
and e'*’® is the general relativistic metric function, k is the total thermal conductivity, Qi, is the 
total neutrino emissivity, is the total specific heat. The first equation is the general relativistic 
definition of photon luminosity and the second equation tells how the photon luminosity varies with 
neutrino emission. 


3.1. Neutrino emission 

In this subsection, we discuss the neutrino emission processes such as the direct Urea, modified 
Urea, nucleon bremsstrahlung, and other processes as summarized in Table Note that the 
neutrino comes out from the entire region of a neutron star before and after its birth. Neutrino 
emission due to PFB will be discussed in Sect. [H 


3.1.1. Direct Urea process 

The most efficient neutrino emission process is the direct Urea process : 


p + e + t'e , 
p + e~ ^ n + Ue . 


(24) 


With charge neutrality, the critical condition (momentum conservation) for the direct Urea process 
can be obtained by neglecting the neutrino momentum, 

PFn = PFp + PFe , 

P%p = PFe + PPti ) (25) 

PFp = . 


^ In the RMF model, we adopt the natural unit c = h = 1, but for the cooling formalism, we follow the convention 
to use the cgs units. 
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Table 6: Neutrino emission processes in the core and crust of neutron stars. Emission rates were 
taken from Ref. |48j . The reduction factor TZ was inserted to account for the pairing effects, n and 
p are quasi-particles {np mixed states). Orders of Qpi and Qpair can be found in Ref. [15] . 



Name 

Process 

Q (erg 

cm ^ s ^ ) 


Direct Urea 

n ^ p + e + Ve 

p + e ^ n 

Qd ~ 

4.0 X 10277^T| 


Modified Urea 

n + n^p + n + e + Ve 
p + n + e^n + n + Vf. 
n+p^p + p + e + De 
p + p + e^n+p + Ve 

y^nn 

Vm 

Vm 

8.1 X 

4.1 X 102l7^T| 



n + n^n + n-\-v-\-v 

/Onn 

7.5 X lO^'^TZT^ 

Core 

Bremsstrahlung 

n+p^n+p-\-v-\-v 

Q7- 

1.5 X 102O7^T| 



p+p^p+p+v+v 

Of- 

7.5 X 10l®7^^| 


Cooper pair decay 

n + n ^ (nn) v + D 
p + p ^ {pp) + v + D 

Qpbf 

- 102l7^^J' 


nn Bremsstrahlung 

n + n^n + n + iy + D 

yonn 

7.5 X 10‘^^TZT^ 

Crust 

eZ Bremsstrahlung 

e + A) —y e A) y D 

Qb- 

lO^'Tgf 

Plasmon decay 

7 —)■ n -t P 

Qpi 



Annihilation 

e~ + ^ V + V 

Qpair 





Figure 5: (Color online) Left : schematic figure for the Urea process. The direct Urea process 
happens near the fermi surface of neutron, proton, and electron. The sum of magnitude of proton 
and electon (or muon) should be equal to or greater than the one of neutron. Right : critical 
fraction of proton to total number of baryons for the direct Urea process. If there exists muon, the 
critical fraction of proton becomes a function of density. 
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Figure 6: (Color online) The left figure shows the proton fraction for the EoS of Skyrme force 
models. In SV EoS, the direct Urea process for electron and muon turns on at lower baryon 
density. On the other hand, SLy4 does not allow both the electron and muon direct Urea process. 
On the right figure, we show the results with RME models. All curves stop when the density reaches 
the central density of the maximum mass of neutron stars. 


If there is no muon in the core of neutron star, the critical fraction of proton for the direct Urea 
process is simply Yp^crit = g- For high dense matter where muon exists and mp, 


Yp^crit — 


1 + ( 1 + ^ 


0.1477. 


(26) 


In general case where both muon and electron exist, the critical fraction depends on baryon number 
density. In Figure we summarize the critical proton fraction for the direct Urea process. The 
dashed line indicates the proton fraction when muon does not exist. The solid line represents the 
critical fraction for the electron direct Urea process. On the other hand, the dotted line shows the 
critical fraction for muon direct Urea process. In Figure we summarize the proton fractions from 
various EoS in comparison with the critical proton fraction for the direct Urea process. In Table 
we summarize the critical densities for the electron and muon direct Urea process. 

The emission rate was firstly calculated by Lattimer et al. |l9] and modified by Yakovlev et al. 
[48] adding superfluidity effects; 


Qd 


= 4.0 X 10 


27 ^n'mp 

mnmp 



1/3 

Tq &npe erg cm“^ s“\ 


(27) 


where m* {m*) is the effective mass of a neutron (proton), Tg = T/IO® K, TZd is the reduction 
factor from superfluidity, and Qnpe is dehned as, 


0 


npe — 



if PFn, PFp, PFe Satisfy the momentum conservation, 
otherwise. 


(28) 
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Table 7: Critical densities {ric in fm“^) for the electron and muon direct Urea process and the 
maximum mass of neutron stars (Mmax) for each model. Numbers in the parentheses {Merit in unit 
of Mq) correspond to the neutron star masses with which the direct Urea processes start to occur. 


Model 

Tic {Merit) 

Yfmax {Mq) 

e-dUrca 

/i-dUrca 

SLy4 

- 

- 

2.07 

SkI4 

0.502 (1.63) 

0.582 (1.83) 

2.19 

SGI 

0.492 (1.72) 

0.616 (2.00) 

2.25 

sv 

0.253 (0.97) 

0.315 (1.30) 

2.44 

TOV min 

0.385 (1.12) 

0.458 (1.37) 

2.05 

LS220 

0.433 (1.31) 

0.527 (1.55) 

2.04 

lU FSU 

0.611 (1.77) 

0.900 (1.94) 

1.94 

DD-ME<5 

0.764 (1.79) 

0.894 (1.89) 

1.96 

SFHo 

- 

- 

2.06 

NLp 

0.340 (1.11) 

0.414 (1.37) 

2.09 

TMA 

0.286 (1.14) 

0.358 (1.43) 

1.99 

NL3 

0.205 (0.85) 

0.255 (1.36) 

2.78 


3.1.2. Modified Urea 

The second important neutrino emission process is, so called, modified Urea proces, 


n + n—)-p + n + e + i^e, p + n + e —>-n + n + z^e) 

n + p^p + p + e~ + He 1 p + p + e~^n + p + Ve. 


(29) 


Compared with the direct Urea process, there is no threshold fraction for protons thus it always 
happens in the core of a neutron star. The emission rate for the modified Urea process is calculated 
by Friman and Maxwell [50] and recalculated by Yakovlev and Levenfish m to add the superfluidity 
effect 


= 8.1 X 10^1 t! anfinTlZ 

V UT-n / V UT-b J \nn J , ^ 

*3 / 1/3 

Q'S = 4.1 X ergcm-3s-> , 

where On and Op are correction terms to describe the momentum transfer dependence of the matrix 
elements in the Born approximation, and fdn and /3p are corrections to the non-Born and to NN 
interaction, which cannot be explained by one pion exchange and Landau theories |51] . We use the 
values of Yakovlev and Levenfish m for an, ap, fin and fip, 


OLn = OLp = 1.76 — 0.63 

fin = fip = 0.68 . 



(31) 


The reduction factors TZ^ and TZ'^ comes from the superfluid nucleons in the core. We used the 
fitting functions provided by Yakovlev and Levenfish m- 
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3.1.3. Nucleon bremsstrahlung 

Emission process important next to the modified Urea process is nucleon bremsstrahlung. The 
difference between the modified Urea process and the nucleon bremsstrahlung is that during the 
process of nucleon bremsstrahlung the flavor of the reaction agent does not change so electron is 
not involved in the process compared with the modified Urea process, 

n + n^n + n + i' + u, 

n+p^n + p + u + u, (32) 

p + p-f-p + p + v + u. 

The emission rate for the nucleon bremsstrahlung was calculated by Friman and Maxwell |50j , and 
Yakovlev and Levenfish m taking into account the superfluidity effect, 

/ \ 4 / \ 1/3 

= 7.5x10^9 ^ ^ TlannPnnKn]Pevgcm-H-\ 

\mnj \no J 

= 1.5 X 10^° f T^anpl3npM^T^'B erg , 

\mnmpj \noJ 

Q'^ = 7.5 X 10^® ^ r|app/3ppAC7^g ergcm“^s“^, (33) 

\mp J \noJ 

where ann, ocpp, otnp are the correction factors from the matrix element in the Born approximation, 
and Pnn, f^pp, flnp are corrections which come from various reasons [5l]. We use the same values in 
Yakovlev and Levenfish 


rT-nn — 0.59, (Trip — 1.06, Oipp — 0.11, 
Pnn = 0.56, Pnp = 0.66, /3pp = 0.7 . 


(34) 


For the number of neutrino flavors, we take AC = 3 including Vt-. TZ^, TZ^, are reduction 
factors and given in Yakovlev and Levenfish m- 


3 . 1 . 4 . Neutrino emission processes in the crust 

In the crust of neutron stars, free protons do not exist and the neutrino emission occurs differ¬ 
ently. In this work, we consider nn bremsstrahlung, eZ bremsstrahlung, plasmon decay, and e^e"*" 
pair annihilation in the crust. 

In the inner crust, free neutrons are available and they are involved in the nn bremstrahlung. 
The emissivity is the same as in the core except for the appearance of filling factor since the heavy 
nuclei in the crust has finite size , 


where / 
given as 


/ ^ \ 4 / \ 1/3 

gr = 7.4 X 10^9 ^ ^ r| ann Pnn A 7 , (1 - /) erg cm -3 s"' 

\mnj \no J 

is the volume fraction of the heavy nuclei in the Wigner-Seitz cell. In the 

mm, 


ann = 1 - -liarctan ( - ) + 777:——77 , 

uj 2{l + u^) 


m^rC 


, (35) 

crust ann is 

(36) 
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Since we assume that there exists a heavy nucleus at the center of the Wigner-Seitz cell, a 
neutrino emission process might happen through the collision between electron and the heavy 
nucleus, 

e + {Z, A) —>■ e + (yZ, A) -\- u -\- u . (37) 

The above neutrino emissivity was calculated by Kaminker et al. and the htting function is 
given as 

QeZ ^ iQh{T,r,p) ergcm-3 3 -^, (38) 

with 

h{T, r, p) =11.204 + 7.304r + 0.2976r - 0.370r2 + 0.188rr - 0.103r2 
+ 0.0547 tV - 6.771ogio(l. + 0.228p/po), 

where r = logj^o^s, r = logj^oPiS) ( 7’12 = p/lO^^ gcm“^), and po = 2.8 x lO^'^gcm"^. 

Though a free electron cannot emit neutrino-antineutrino pair due to the lack of energy- 
momentum conservation, electrons with interacting medium can emit neutrinos via e ^ e + ly + u. 
The process can be written as 

j ^ u + u , (40) 

where 7 stands for a plamon. This process is highly efficient at high temperature and medium 
densities in the crust of neutron stars |48j . Simplihed ht is given in Yakovlev et al. 


(39) 


Qpi — 


Qc 


967r^a V T, 


(16.23/6 + 4.604/;-6)e-7r ^ d 


where 


G% /mpC\9 




h \ h 


1.023 X 10 ^^ergcm"3s“S ^(7^ = 0.9248, Tr = 


nieC 

ks 


(41) 


5.93 X lO^K, 


and fp is the electron plasma parameter defined as 




_ hupe 4:Tre‘^ne/ml 


keT keT 

e“e+ annihilation might happen in the crust of neutron stars, 

e~ + e~^ ^ V + V. 

The formula for e^e"*" annihilation is given in Eq. (22) in Ref. 

3.2. Heat capacity 

The specific heat in neutron star is given by the sum of its constituents. 


(42) 


C, = Y,Ci 


(43) 


where i denotes the type of particles. In general, Ci is given by 

m*pFik‘^T 


Ci = 


3h^ 


(44) 
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where m* is the effective mass at the Fermi momentum ppi. For neutron star core, we have the 
contribution from neutrons, protons, electrons, and muons. 


= Cn + Cp + Ce + C^. (45) 

Instead of using the simple formula for the specific heat of electrons [53] , we use analytic formula 
for leptons in Lattimer-Swesty m as discussed below. The pressure, energy density, and entropy 
density contribution from leptons are simply functions of chemical potential. For example, the 
energy density is given by 


^ ^ / W 


3 r 


1 + ^ 


\hc 

where gi is the degeneracy factor for lepton and pi is a chemical potential of lepton. 


1^1 = r 


r = 


■\/q^ +1"^ +t 


ll/2 


(47) 


(48) 


where t = {hc)^ni/gi and q = {7rT)‘^/3 — mfc^f2. The electron specific heat can be obtained by 
evaluating dSi/dT. 

For neutron star crust, there is no contribution from protons but there are contributions from 
lattice ion, 

= C, + Cn + Cion. (49) 


Depending on the temperature, the ions are in the gas, liquid, or solid phase. The F factor which 
is the ratio between Coulomb energy and thermal energy, tells us in which phase the ionic lattice 
is as discussed in Ref. [54] . 


ksTvc 


(50) 


where rc is the radius of Wigner-Seitz cell {^Trr^ni = 1). The contribution from the lattice ion 
becomes 


c. 


3nikBfD{T/9D) if 


where the internal energy of the Coulomb liquid (F < 
given by (SHE] 


u 


riiksT 


= F^/^ 


A^ 


+ 


V A 2 + F 1 + F 


We use values of ^ 1 , 2,3 and i?i, 2 , 3,4 as in the Ref. |6]. 


F < 1 (Gas Phase) 

1 < F < 178 (Liquid Phase) (51) 

F > 178 (Solid Phase) 

Fm, where F^ is the melting constant) is 


^ B2 + r ^ S 4 + r4 • 


(52) 


0£) is a Debye temperature 6 d = 5.6 x 


^ A simplest formula for the specific heat of electrons |53] is given as 


De 


m*pF^k%T 

3?P 


5.67 X 10^® 



2/3 




erg 

cm^Jf ' 


(46) 
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Table 8; The best fit parameters for ^’s and -B’s in Ref. [6]. 



^2 

CO 

B2 

B3 

Bi 

-0.9070 

0.62954 

-\/3/2 - Al/^/Z2 0.00456 

211.6 

-0.0001 

0.00462 


• fm^/^)^/^Rr. fo is defined by 


fnix) = 3x^ 



dt 

(e* - 1)2 ’ 


and we use simple analytic fitting function [56] as 


fnix) 


-0.2172+0.05014x+23.0588x^ 

1.+23.0588x2 
1 - 0.05X-2 

\ 


if x<0.15, 
if 0.15 <x <4, 
if 4 < X . 


( 53 ) 


(54) 


For the smooth transition from liquid to solid phase, we use linear combination when T is between 
T = 178 and T = 220 [57], Cion = C'l(T)^^^^ + Cs{T)^-^^^ where Cl{T){Cs{T)) is the specific 
heat of liquid (solid) phase. 

The main contribution to the heat capacity comes from the core since it has large volume and 
mass fraction of a neutron star. If there exists nuclear superfluidity, Cn and Cp become highly 
suppressed as temperature decreases to much below the critical temperature [53]. In this case, the 
electrons (or muons) are the main contributor to the heat capacity of a neutron star. Thus the heat 
capacity is 20 times lower than the one without superfluidity since the electrons are not suppressed 

m- 


3.3. Thermal conductivity 

The thermal conductivity arise from the collision phenomenon between particles for given den¬ 
sity and temperature. The kinetic equation for thermal flux density is given by [591160| 


qi = 


1 

47r3 


dki {a - Fi 


-KiVT, 


(55) 


where kj, Vj, Ci, fii are the wave-vector, velocity, energy, and chemical potential of particle i, 
respectively. Ki is the thermal conductivity and T is the temperature. Fi is the distribution 
function [60] 

Fi = (56) 

where fFD,i is the Fermi-Dirac distribution function and is the deviation from the equilibrium 
distribution. This can be obtained from Boltzmann equation, 


= -nici - m)- 


■ VT 


(57) 


where r* is the relaxation time (inverse of the collision frequency, I/tj = t'j) and should be calculated 
numerically using collision integral. Then the thermal conductivity has the compact form |60j . 


Ki = 


ir'^k'^TniTi 

3m* 


(58) 
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where n* is the number density and m* is the effective mass. 

In the core, the thermal conductivity consists of neutron, electron and muon contribution. 




( 59 ) 


Since a neutron can collide with both neutron and proton, the relaxation time is given by 


Tn 


1 

^nn “ 1 “ ^np 


(60) 


The electrons also collide with proton, electron, and muon. Since there are two different lepton 
species (electron and muon) in the core, the collision relaxation time and frequency relation is 
different from the single lepton (electron) matter |61j : 


1 = VeTe + 1^'epTp, 

1 = , 


(61) 


where the prime indicates the particle after collision. In Vf,, we have to consider the collision 
frequencies from ee, e/x, and ep, 

Z^g — Z/gg T T l^ep • (1*^) 


In the core, we use the results from Ref. [60] for neutron thermal conductivity and Ref. |61j for 
electron and muon thermal conductivity. Both papers consider the superfluid nucleon state so it 
gives reduction factor for a given temperature. In the crust, electrons are the main thermal con¬ 
ductivity factor and an electron collides with an electron or an ion. Thus, the thermal conductivity 
is given by 


He 


TT^r/c^ng 
3m* Ue 


Z^e — Z^ee T ^ei ; 


(63) 


where m*c = ymgC^ Ppe- result from Ref. |62| for electron thermal conductivity 

in the crust of neutron star. Since the critical temperature for electron superconductivity is much 
lower than the internal temperature of neutron stars, it is not necessary to consider the superfluidity 
reduction. 

As in the case of heat capacity, once the temperature drops down the critical temperature of 
nuclear superfluidity, the thermal conductivity caused by collisions between superfluid baryons or 
between electron, muon and superfluid baryon experience reductions. Thus the electron thermal 
conductivity dominates both in the core and crust. 


4. Results from the Standard Cooling 
4-1- Standard cooling and direct Urea 

Simulation shows that the mass of most of the compact remnants after the supernova explosion 
is below 1.5Mq [63], and similar distribution of mass range is obtained from the precisely measured 
masses in NS-NS binary systems. Therefore, to a good accuracy, we can assume that most of the 
sources given in Table [^ have masses less than 1.6Mq or 1.7Mq. 

Figures [7] and [^ show the cooling curves for Skyrme models. SLy4 model shows similar behaivor 
regardless of the mass of neutron stars, which is mainly due to the absence of the direct Urea process. 
Were it not for the direct Urea, modified Urea is driving cooling mechanism in the standard cooling 
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scenario. We can see that the modified Urea is good at reproducing the data for young (below 10^ 
yrs) and old (above 10® yrs) stars, but completely misses the middle age (10^ ~ 10® yrs) data. This 
may imply that actual cooling will go through slow-quick-slow stages of neutrino emission process. 
For SkI4 and SGI models, temperature drops abruptly for masses 1.7Mq and l.SM©, respectively. 
This abrupt decrease of temperature is the signal for the ignition of the direct Urea. In fact, the 
direct Urea is the fastest neutrino emission process ever known, so once it is turned on, regardless 
of the existence of PBF or exotic states, shape of the cooling curve is predominantly controlled by 
the direct Urea. The result shows that the direct Urea is too fast that it fails to pass through any 
observation data. On the other hand, if cooling is driven by the modified Urea, SLy4, SkI4 and 
SGI show very similar thermal evolution trajectory. In the SV model, the direct Urea occurs in 
all the neutron stars, so the model cannot explain the temperature data at all. For TOV min and 
LS220 models, the modified Urea is the main mechanism for low mass stars, but the direct Urea 
starts to occur also in the low mass stars. Assuming that most of the mass of the measured star in 
Figurej^is in the range of I.OMq ~ 1.6Mq, TOV min can hardly explain the observed temperature 
profile. It is striking that though TOV min model shows similar quality of mass-radius relation to 
SLy4, SkI4 and SGI (Figure]^, they predict quite different thermal evolution scenario. Combining 
the empirical data from both mass-radius relation and temperature, we can reduce the space for 
nuclear models which are applicable to the investigation of superdense nuclear matter. For this 
reason, we remove SV and TOV min models from the consideration hereafter. 

Figures and [T^ present the cooling curves with RMF models. Now we can easily see that 
the direct Urea is not working in the SFHo model, and it is activated only in large mass stars in 
lU-FSU and DD-ME5 models. We note that mass-radius relations for SFHo, lU-FSU and DD-ME5 
models are similar to those of SLy4, SkI4 and SGI models. The similarity is kept for the cooling 
curves, which may evidence a strong correlation between bulk properties of neutron stars and their 
thermal evolution. Three hard EoS models, NLp, TMA and NL3 cannot reproduce the observation 
data at all, so we exclude them in the coming analyses. 

A feature worthy of notice is that the cooling curve due to the modified Urea is sensitive to the 
fraction of protons in the core of neutron stars. This can be easily seen by comparing the cooling 
curves of SLy4, SkI4 and SGI, which are similar to each other, to those of TOV min and LS220. 
The latter models show fast cooling at young (less than 10^ yrs) and middle (10^ ~ 10® yrs) ages, 
but the temperature drop slows down significantly at old ages (above 10® yrs). One can observe 
a similar patten from the results of RME. This distinctive behavior in the modified Urea cooling 
can be partially understood from the proton fraction in Eigure With more protons (TOV min 
and LS220 in Skyrme force models and NLp and TMA in the RME models), we have cooler stars 
in the young and middle ages. 


4-2. Radius and symmetry energy properties 

The radii of neutron stars have a close relation with the pressure around nuclear saturation 
density [65], R oc The energy per baryon and the pressure around saturation density can be 

expanded as 


E{n,d) 

P 


-B + 


V 3 no 



L 

3 no 


52 



(64) 

(65) 
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Figure 7: (Color online) Surface temperature (left column) and photon luminosity (right column) 
vs age without superfluidity effects in non relativistic SLy4, SkI4, and SGI models. The symbol 
indicates the effective temperature of Gas A neutron star. Each curve in the plot corresponds to 
different neutron star mass in the range of I.OMq to 2.0Mq. In case of SLy4, the direct Urea is 
not turned for any mass of neutron stars. We use the (GPE) Tg — Tb relation in |64) . 
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Figure 8: (Color online) Same plot as in Figure]^ for SV, TOV min, and LS 220 models. SV has 
high proton fraction even in the low mass neutron stars so it shows the fast cooling processes. In 
case of TOV min, the direct Urea process happens even in 1.2M0. GPE Tg — T}, relation was used. 
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Figure 9: (Color online) Same plot as in Figure]^ for relativistic mean field models; lU-FSU, DD- 
ME<5, and SFHo. For RMF models which have maximum mass less than 2.0Mq, the cooling curve 
starts from I.OMq and end up with 1.9Mq. For others, we draw up to 2.OM0. Critical neutron 
star mass for the direct Urea process in lU-FUS model is \.77Mq but the significant effects are 
only for 1.9Mq neutron star. CPE Ts — Tf, relation was used. 
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Figure 10: (Color online) Same plot as in Figure for NLp, TMA, and NL3. The direct Urea 
process for NL3 happens even less than I.OMq so that all curves show fast cooling. GPE Tg — Tb 
relation was used. 
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where 6 = - 
energy (L), 


—. This relation also gives a rough relation with the density derivative of symmetry 


R oc 


( 66 ) 


The proton fraction is determined from the ground state energy of nuclear matter. The symmetry 
energy, which is a qualitative estimate of attraction between neutron and proton, has the relation 
with proton fraction. The algebraic relation S{n) = Sy ^ indicates that the greater L, the 
greater proton fraction. This implies that the direct Urea process is connected with the radius of 
neutron stars even though the proton fraction at high density depends on the higher order terms 
of symmetry energy, anc 


neutron stars. Figure 
stars’ radii would tel 


11 


thus L is still a good indicator of the direct Urea process in the core of 
shows that is indeed model independent. Roughly speaking, neutron 
us the amount of proton fraction so, we can guess that the direct Urea 
process would be turned on or off. Thus, too large radius {RiaMq > 14 km, or L > 90 MeV) is 
not favored since the direct Urea process happens even for the small mass of neutron stars. This 
is also consistent with Steiner et al. [3] in which they estimated the area of neutron stars’ masses 
and radii using X-ray burst data. 

Lattimer and Lim [66] summarized symmetry energy properties {Sy, L) both with experimental 
results and theoretical calculations.The analysis from the nuclear mass fits, neutron skins, heavy-ion 
collision, giant dipole resonances and dipole polarizabilities gives an overlapped region. Considering 
the theoretical calculation of pure neutron matter and astrophysical observations of neutron stars, 
the allowed ranges of symmetry energy (Sy) and its density gradient (L) are 29.0 MeV< Sy < 32.7 
MeV and 40.5 MeV< L < 61.9 MeV. Our result for neutron star cooling indicates that L < 85 MeV 
so that the direct Urea process should not be turned on in the low mass neutron stars (M < 1.2Mq). 
This is consistent with Lattimer and Lim’s conclusion. 
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Figure 11: (Color online) RiaMqI various nuclear models. RiaMq /is independent of 
models. 
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Figure 12: (Color online) Band plot of from light elements envelope and heavy elements 
envelope. Each band has the mass range between 1.2 Mq and 2.0 Mq. SLy4 cannot explain some 
of data even with the union calculation of light and heavy elements since the direct Urea process 
is not activated even in the maximum mass of a neutron star. On the other hand, SkI4 EoS can 
cover all the observed data. 


4-3. Effect of envelope elements 

The surface temperature highly depends on the chemical abundance of light elements [ISIET]. 
Figure [T^ shows the band plot of both with the light elements and heavy elements. The bands 
in each plot indicate neutron star masses in the range of 1.2Mq ~ 2.0Mq. At early age, the top 
curve represents the most massive star and later time, the curve from massive stars is the one with 
the lowest temperature. 

For the real cooling process, one has to take into account the cooling from both light and heavy 
elements. The chemical evolution from the light elements to the heavy elements or pulsar injection 
of light elements into the magneto sphere |16j indicate that the real cooling curves may start from 
the band with light elements and moves towards the band with heavy elements as the neutron star 
evolves. The mass of light elements is dehned as 

AM{t) = AM{t = . (67) 
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Figure 13: (Color online) Light element decay and neutron star cooling curve. Left panel shows 
1.4Mq neutron star cooling path with SLy4 EoS. The initial mass of light elements is assumed to 
be 10 “'^Mq and the decay starts at t = 10^ years. The right figure shows the same evolution path 
but with two different masses of neutron stars (1.4Mq, 1.7Mq) using SkI4. 


where Td is the decay time. In the accreted envelope, the surface temperature can be fitted as a 
function of the mass fraction of light elements to the total mass of neutron stars 167]. If the surface 
is made of pure irons, 

T^S(>,Fe = gum?-"'’ + iC/m"'] , ( 68 ) 

where ( = T^g - ( 7 rb 9^/^)^/^/103 and 514 = (l - We define Tb (Tbg = 

Tb/10®K) is the temperature where the energy density is 10^*^ g cm“^. On the other hand, if the 
surface is fully accreted, i.e. there are only hydrogen, 

7;W = 514(18.1 rb9)2-42. (69) 

For partially accreted envelope, with the definition of rj = g^^AM/M, we have Tg — Tb relation. 


Tg = 


_l_ 7"4 1 1/^ 

“-‘efr6,Fe -^efr,a 

Cl “h 1 


(70) 


where a = [1.2 + (5.3 x 10“®/r/)°'^®] Tj^g^. Using Tg — Tb relation in Ref. [HTj, we can find the surface 
temperature for given mass of light elements on the surface. 

In Figure [T^ we show the cooling paths obtained by taking into account the light element decay. 
Depending on the amount of light elements and the decay time scale (r^), the actual neutron star 
cooling curve will be located between the light and heavy elements bands. In this figure, the initial 
mass of light elements is assumed to be 10 “'^Mq and the light elements start to decay when t = 10 ^ 
years after the birth of neutron stars. In case of SLy4, the cooling curve is almost identical for all 
mass of neutron stars, thus a typical mass 1.4Mq is chosen to see the evolution path. For SkI4, 
1.4Mq and 1 . 7 M 0 cooling paths are shown to compare the effects of elements decay. This implies 
that elements composition, decay history, and the direct Urea process can be used selectively to 
explain all the observed data. Note that the EoS which doesn’t allow the direct Urea process (i.e. 
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Figure 14: (Color online) The direct Urea Process effects on LS200 and lU-FSU. For both models, 
most of the observations are in the very narrow mass range (AM = 0.01 Mq). GPE Tg — Ti, relation 
was used. 


SLy4 for example) cannot explain mid-old low temperature neutron stars without other fast cooling 
mechanism such as cooper pair emission or Bose condensation. The paths with light element decays 
also show that the rapid temperature drop can occur during the decay process of light elements. 

5. Results with Superfluidity 

It is believed that ^S'o neutron superfluidity exits in the inner crust of neutron stars, and 
proton and ^P 2 neutron superfluidity state appear in the core of neutron star. Since there is no 
free proton in the crust of neutron star, the superfluidity involving proton exists only in the core of 
neutron star. Beta equilibrium condition allows the fraction of proton only up to 0.2 so the effect 
of the proton superfluidity itself is not as strong as that of the neutron superfluidity. However, the 
proton superfluidity delays the surface temperate drop mainly due to the reduction factors in the 
neutrino emissivity, heat capacity, and thermal conductivity. Once the local temperature drops 
below the critical temperature for superfluidity, the properties of nuclear matter change drastically. 
The neutrino emissivity, heat capacity, and thermal conductivity are reduced significantly and 
behave like ~ exp(—A/T) where A is a paring gap energy. The numerical reduction factors for 
each physical quantities are given by Yakovlev et al. [Ml EU]. As discussed in Page and Applegate 
|68j and Yakovlev et al. [53] . the direct Urea process turns on in the really narrow mass range. 
As shown in Figures - [Tol the direct Urea process imposes huge effects on the cooling curve. 
If the mass of neutron star is slightly greater than the critical mass for the direct Urea process, 
(M > Md + O.OIMq), then the direct Urea effect is clear 

In Figure [T^ the effects of the direct Urea process on neutron star cooling are summarized with 
two EoS LS220 and lU-ESU. The results show that most of the observed data can be explained by 
two curves 1.31Mq and 1.32Mq for LS220 and 1.78Mq and 1.79Mq for lU-ESU. Thus, if there’s 


^Numerically, there needs more grid points in the core of neutron star to treat the direct Urea process properly. 
In this case, we used 16 times more grid points to see the split of the curves between I.SOMq and l.dOM©. 
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no superfluidity effects, the mid age cold neutron stars (e.g. Vela, Geminga) should have the mass 
in the really narrow range. In the statistical point of view, it is unacceptable that almost all the 
observations are in the 0.01Mq/Mq ~ 1% since the mass distribution of the observed neutron stars 
has the broad range [TO]. This problem can be managed if we employ the pairing effects. The 
neutrino emissivity formula for PFB is given by |69l 170] 


= 3.51 X 10^^ 



PFi \ 

rriiC J 



erg 
cms ’ 


(71) 


where i represents type of nucleons {i = n,p) and j stands for singlet (j = s) or triplet {j = t, mj = 
0) pairing. Fg and Ft are given in Ref. [69j 



dx 

(1 + ' 



dx 

(1 + ’ 


(72) 


where x is an integration variable, y = Aj(T)/T, z = and f dfi is the solid angle 

integration. The fitting functions for Fg and Ft are also given in Ref. |69) . aij’s are given by Ref. 

izni 




(73) 


0,n,t — , Up^t — p , 

where Cv,n = 1, CA,n = 9A, Cv,p = 4sin^0rv ~ Ij Ca,p = —QA with qa ~ 1-2 and slv? 9y/ zz 0.23. 
Since the core of neutron stars is mostly composed of neutrons {ppn ^ PFp) and the magnitude 
of On,t and Op^s are comparable, the triplet PBF is the main neutrino emission agent once the 
superfluidity occurs. (Note that ^P 2 neutron and proton pairing are expected in the core of 
neutron stars.) When the temperature drops below the critical temperature, the modified Urea and 
bremsstrahlung neutrino emission processes are highly suppressed and PBF process overwhelmes 
the other neutrino emission m- 

Aside from the previous calculations of nuclear superfluidity, we introduce the phenomenological 
pairing gap formula to see the effect of gap size and the density range. 


Tcikf) 


-M-ikf- koTYk2 - kff- if ko<kf< fe; 
0 if otherwise. 


(74) 


where is the maximum critical temperature for superfluidity for given kf (the fermi wave 

number for total bayron number density), /cq (^ 2 ) is the starting (ending) wave number for given 
type of pairing. M is the normalization factor for the critical temperature. 


Af = 




1 


(75) 


Most of the neutron gap calculations can be fitted with ac = 2 and f3c = 2. Since the exact 
proton fraction in dense nuclear matter is not certain, the above formula is a good approximation 
to check the effects of the pairing on the specihe heat and neutrino emission. 
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Figure 15: (Color online) Critical temperature for different type of pairing, kf is the wave number 
of Fermi-momentum (in unit of fm~^) for the total baryon number density. The critical temperature 
is then evaluated for beta-stable nuclear matter. The core-crust boundary is obtained assuming the 
phase transition happens at p ~ l/2po = 0.08 fm“^. Each type of pairing can be obtained using 
different methods. ^Sq neutron : CBF - Correlated Basis Function m, PP - Polarization Potential 
|73| . BHF - Brueckner Hartree Fock m. RG - Renormalization Group j75j. ^Sq proton : DBHF - 
Dirac Brueckner Hartree Fock [76], PCT - Parameterized Critical Temperature |77|, BHF |74] . ^P 2 
neutron : BHF (TS), PCT [77], OPEC (BCS) - One Pion Exchange Gaussian with generalized BCS 
m- Note that proton superfluidity is not operative in the crust of neutron stars since free gas 
of proton is not available in the crust. The right bottom figure shows the critical temperatures for 
each type of pairing. All curves were obtained from the simple Eq. (74) with ac = 2 and (3c = 2 for 
and ^P 2 neutron pairing and with Uc = Q and (Sc = 1.2 for proton pairing. 
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Figure 16: (Color online) The effect of superfluidity for SLy4. 


In Figure 15 


the critical temperatures are summarized for given kf in beta-stable nuclear 
matter. Thus for ^P 2 neutron pairing (left bottom), kf may be different from the original paper 
since we use kf for the total baryon number density not for the pure neutron matter. If the gap 
calculation is done in the pure neutron matter {kpn), the proton fraction is given by APR EoS [71| 
to recover kf oi the total baryon number density for this hgure. 

In each gap calculation, the critical temperature strongly depend on the methodology. Consid¬ 
ering this fact, we use the phenomenological critical temperature formula and see the cooling curve 
how it depends on it. For and ^P 2 neutron pairing, Uc = 2 and Pc = 2 are suitable to represent 
the critical temperatures. For proton pairing, we adopt Uc = Q and Pc = 1-2 to mimic the 
behavior of the critical temperature in this example. 

To see the effect of nuclear pairing, we choose SLy4 and SkI4 to simulate neutron star cooling, 
because both SLy4 and SkI4 satisfy the mass-radius constraint zone [3] and SLy4 does not turn on 
the direct Urea process, while SkI4 does if the mass of neutron star is greater than 1.63 Mq. The 
cooling curves for the case of SLy4 in Figure 16 show that the early start of ^P 2 pairing gives the 
narrow band of cooling curves. The early start means ko is close to the crust and core boundary. 
Thus, the core of neutron stars has the superfluidity regardless of its mass. Therefore, small ko 
gives sharp drop of temperature in young age of neutron stars. In this case, the observational data 
for old age neutron star cannot be explained with cooling simulation. On the other hand, if the ^P 2 
neutron pairing appears in some range of density, the low mass neutron stars do not show sharp 
temperature drops since ^P 2 pairing is not available in the core of neutron stars. It seems that 
if the EoS does not allow the direct Urea process, the combined area of both from the light and 
heavy element envelope cannot explain the mid age low temperature neutron stars or old age high 
temperature neutron stars at the same time. This might indicate that pion or kaon condensation 
might exist in the core of neutron stars if the massive neutron stars does not have enough proton 
fractions. As summarized in Eigure [T^ since SkI4 allows the direct Urea process, the bands from 
cooling curves can explain every observational data without pairing effects. If we assume pairing 
effects, however, reduction factors for heat capacity and neutrino emissivity play a role to restrict 
ko and k 2 for ^P 2 neutron pairing. The critical temperatures for SkI4 are summarized in Table 
For each case, we use Oc = 2 and Pc = 2 for ^S'o and ^P 2 for neutron, and Oc = 8 and Pc = 2 for 
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Figure 17; (Color online) Superfluidity effects with SkI4 EoS. Each case has the different critical 
temperature so different dlnT/dlnf. The error bars in the right bottom plot denote the analysis of 
ACIS-S (Graded Mode) in Ref. [80]. The solid lines in each plot show the cooling curve of 1.5M© 
neutron star with AM = 5 x 10“^^M©. 
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Table 9: The critical temperature parameters for each case. 




Case I 

Case II 

Case III 

^P 2 n 

jnmax 

ko (fm-^) 
k 2 (fm“^) 

6.80 X 10^ 
0.99 

2.8 

6.47 X 10® 

1.5 

3.0 

5.95 X 10® 

1.8 

2.4 

^Sop 

jnmax 

ko (fm-^) 
k 2 (fm“^) 

6.44 X 10® 

0.1 

2.5 

7.45 X 10® 

0.1 

2.5 

1.0 X 10® 

1.8 

2.4 

^S'o n 

j^max 

ko (fm"^) 
k 2 (fm“^) 

3.2 X 10® 

0.0 

1.3 

5.0 X 10® 

0.0 

1.3 

1.0 X 10® 

0.0 

1.3 


for proton critical temperatures. The exact critical temperature might be obtained if we have 
the exact temperature drops of the Cas A neutron star for the next decade. 

6. Conclusion 

In succession to our prior works on the mass-radius relation of neutron stars [laiia], in this work, 
we investigated the consistency of nuclear models with observation of neutron star temperatures. 
First, model selection was performed by constraining the model prediction of maximum mass of 
neutron stars larger than 2Mq. As a result, we picked up 6 models among the non-relativistic 
Skyrme force models, and 6 models among the RMF models. In the second step, with the 12 
models selected, we calculated the cooling curves with only standrad cooling mechanisms and the 
direct Urea processes. The result shows clear dependence on the EoS of a model. If the direct Urea 
is not operating, the modified Urea is the most efficient way to emit neutrinos. We found that the 
modihed Urea reproduces the observation data well for the ages less than 10^ years or more than 
10^ years. However, no model could explain the data in the age of 10^ ~ 10^ years. The modified 
Urea always gives temeperature much higher than the observed ones. Once the direct Urea is 
working inside of neutron stars, the star cools down so fast that the cooling curves never touch 
the observation data. Thus, nuclear models which have large symmetry energy gradient {L > 85 
MeV) cannot demonstrate any of astronomical data for neutron stars since the direct Urea process 
is turned on even for low mass neutron stars. As a result, we could sort out the nuclear models 
that satisfy both mass-radius relations and temperature data. 

Surface temperature, which is a directly measured quantity, heavily depends on the composi¬ 
tion of elements in the envelope. Models with moderate symmetry energy gradient can explain 
observational data with standard cooling processes. In the middle age, high temperature neutron 
star indicates that the surface of neutron star has some composition with light elements. With the 
decay of light elements to heavier ones, neutron star cooling curve can be any of one in the both 
light and heavy elements bands. Thus if EoS bands from light and heavy elements contain all of 
the observed data, it can be a good model for the neutron star cooling simulations. 

We have explored the effect of nuclear pairing by combining the PBF to models that are qualified 
with the standard cooling mechanisms. We could reproduce the middle age low temperature data 
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which are missed in the standard cooling mechanism, and at the same time could satisfy the Cas 
A temperature profile. However, looking into the detail, we find that the cooling curves with PBF 
are sensitive to the physical inputs such as pairing gaps, critical temperature for PBF, neutron 
star mass, and etc. Moreover, the temperatures predicted from the curves that are consistent with 
the middle age and Cas A data are too low for the stars with ages more than 10® years, so the 
curves with PBF scarcely touch the old age data. Measurement of temperature change of Cas A 
in the next decade will shed some light on resolving these problems and we expect to reduce the 
uncertainties in the underlying physics. 

In conclusion, we could confirm that the existing mass-radius relation and thermal evolution 
history of neutron stars provide critical test grounds with which we can narrow down the space 
of nuclear models which are good candidates for the application to dense nuclear matter. In this 
work, we concentrated on nuclear models and uncertainties due to envelope elements and nuclear 
pairing. In a recent work m, we demonstrated that hyperons can exist in the core of neutron 
stars which satisfy the empirical mass-radius relations. Urea process accompanied with hyperons, 
though not so efficient as the direct Urea, is another source of fast cooling mechanism. Neutron 
star cooling with hyperons will be coming up in the near future. 
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Appendix A. Numerical Solution for Neutron Star Cooling 

The analysis for numerical solution of the diffusion equations can be found in the appendix of 
D. Page’s thesis [57j. Here we adopt the same notation with Page’s thesis in which tri-diagonal 
scheme is used to solve the coupled diffusion equations. To solve the diffusion equations, Lr (T) 
is defined on the even (odd) grid. We present the numerical solutions for the neutron star cooling 
using penta-diagonal and tri-diagonal schemes and compare the convergency between them. 

Appendix A.l. Penta-diagonal scheme 

The penta-diagonal scheme appeared in ref. |56j . The first diffusion equation is 

C = -PS^ (A.l) 


®In penta-diagonal scheme Lr (T) is defined in even (odd) grid and intermediate time step, and in tri-diagonal 
scheme Lr and T are defined in all grid points. 
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where 

£ = e2'*’«L^, r = e^«r, P = k, S = {Airr^fe'^on. (A.2) 

Here a(r) is the enclosed baryon number in radius r and it is related by 


— = 4:TTr‘^n\l 1 — 
dr 


2Gm 


rc‘^ 


(A.3) 


where n is a baryon number density. The second diffusion equation becomes 


dT 

dt 


-Q^-R 


(A.4) 


where 

(3 = ^, = (A.5) 

Lyy 

As usual, we define C on the even grids (£o, ^ 2 , • • •, T 2 /+ 2 ) and T on the odd grids (71, Ts,..., T 21 + 1 ) 
since we can use the boundary condition at the center {Cq = 0). The main feature of the penta- 
diagonal scheme is that the luminosity and temperature profile can be obtained even in the mid 
time interval. Thus the unknowns for numerical equations are 


/^n+1/2 pn+l /T'n+1/2 ^n+1 ^n+1/2 ^n+1 q-n+1/2 .^.^+1 ^n+1/2 .y.n+l\ 

V -^0 ’•: 'I )-^2 >''^2 > • • • ) ' 2/+1 ’ ' 2 /+ l >''^ 2/+2 ’' 27 + 2 /’ 


where the subscript means spatial dimension and superscript represents time step. The total 
number of unknown are 4/ + 6. Using these mid time interval variables and Henyey method, we 
write the thermal evolution of a neutron star time index from n to n + 1 as 


=nui - A* 


^n+l/2 dC 

da 


n+1/2 

2i+l 


R. 


,n+l/2 

2i+l 


(A.6) 


where 


^n+l/2 _ pn+1/2 Qn+1/2 d7~ 

^2i - -^2i ^2i ^ 


n+1/2 


2i 


(A.7) 


The thermal evolution of time index from n + 1/2 to n + 1 is given by 


Tjs;=- Ut 


dC, 


n+1 

2i+l 


, pn+1 

+ 7^2 *+1 


(A.8) 


where 


^n+l _ _pn+lQn+l((// 

^2^ - ^2i ^2i ^ 


n+1 


2i 


(A.9) 


In general the thermal conductivity k (= P) is a function of temperature T, and the temperature is 
defined only on the odd grid in the numerical scheme, thus we use the average = 5 (T 2 j_i+i-’ 2 i+i)- 
From the above equations and P 2 i, we have equations to solve 


p _ p f ^n+l/2 .pn+1/2 .y-n+1 ^n+l/2\ 

7*4i+l — ’ '2j+l > '2i+l >''^2j+2 > 


_ pn+ 1/2 pjn+ 1/2 

~ -^ 21+1 + ^2i+l 


'2i+l > '^2j+2 
^n+1/2 ^n+1/2 

-^ 2 i +2 ~ ^2i 

da2i + da2i+i 


+ 


q-^+l _ T-n 

^ 2 i+l '2i 

At 


n 

2i+l 


= 0 , 


(A.IO) 
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(A.ll) 


F. 


4i+2 


= F, 


4i+2l 


j rn+^ ,-T-n+l/2 ^n+4 r'ri+4 'i 
\^2i ’ '2i+l ’ '2j+l’''-'2i+2^ 


_ p«+l 
~ -^2i+l 


n+1 

2j+l 


+ Q 

F4*_1 = F4._i(r2r|/' 


rn+l _ rn+l 
‘~'2i+2 ‘~'2i 

da2i + da2i+i 

^Ti-\~4j2 r-r-n-\~4l2\ 
i^2i ) '2i+l ) 


+ 


-T-n+l _ q-‘ 
'2i+l 


■n+l/2 


2i+l 


At/2 


= 0 , 


_ ^n+ 1/2 


2i 


+ 


3 n+l /2 

2j-l 


'2i+l 
3 n+l /2 
2i+l 


•^2* 


/-y—Ti + 1/2 71+1/2 

^ 2 i+l ^ 22-1 

da2i + da2i+i 


77 _ 771 /T-n+l rn+l T-n+lx 

^4i - -t'4ii/2i-H+i > / 2 i+lJ 

pn+1 I pn+1 
-'^ 21-1 “T -'^ 2 i+l 


= + 


Szi 


/ T - n+l _ / T - n+l 

^27+1 ^ 2 i-l 

da2i -|- da 2 j+i 


= 0 . 


= 0 . 


(A.12) 


(A.13) 


These four types of equations can be solved by the multi-dimensional Newton-Raphson method. 
The variations of SC’s and 5F's can be obtained by solving linearized Newton-Raphson method, 


F4*+i + A4.+i<5£^+^/' + Cii+idT^-i^ + Dii+i5T:^Xl + ^4i+i<5T2i/" = 0, 

F4*+2 + Aii+2dCit^ + B^i^2drXt+T + C'4*+257^/J + Eu+2dLltX2 = 0, 
F4*-i + + C4i-1(54/^/' + = 0 , 

F4* + + E^bl^^Xl = 0 , 


where 
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(A.14) 

(A.15) 

(A.16) 

(A.17) 


(A.18) 

(A.19) 

(A.20) 

(A.21) 

(A.22) 

(A.23) 

(A.24) 

(A.25) 
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(A.26) 
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(A.27) 


(A.28) 


(A.29) 

(A.30) 


(A.31) 


In the above equations to solve Fi,..., A 4 /+ 2 , we have unknowns , y^+ 1 / 2 ^ , ..., 

^ 2 i> 2 ^’ and F 2 t +2 ■ There are 41 + 6 unknowns quantities with only 47 + 2 equations, but one 
can obtain the solutions with four additional boundary conditions. With the Tg — Tb relation 
IM1E3EI], we can have 

F 21+2 = e^'^^^+^{47rR^aB)T^ = e^'^^^+^47rR^aB)f{T2i+i). (A.32) 


^Lr{r = 0) = 0 reduces two unknowns 
L 21 + 1 ) reduces and ^ function of 


= 0), and the uniform luminosity approximation (I/ 2/+2 
and respectively. 
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We have matrix equations to solve 


(Cl Di El 0 
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0 ^4 0 (74 
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\ ( \ 

< 577+1 
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A 41-1 
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Ail 
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AI 


D 


4/+1 


^ 47 + 2 / 


dC 


n+ 1/2 




C41-1 0 -E47-1 0 

0 C41 0 

^47+1 0 (747+1 

0 ^47+2 -B 47+2 

(A.33) 

This penta-diagonal linear equation can be solved hy L — U decomposition or Gaussian elimi¬ 
nation method. 


(5£^+i 


F 4 
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5Clt^ 
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F 41+1 
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E 2 

E, 


Appendix A.2. Tri-diagonal scheme with every grid point 

Another method to solve the diffusion equation is to use every grid point {Cq, Ci, • • •, Cn, 
7i, Ti, • • •, Ttv)- In this case, we mix forward and backward numerical differentiation to make 
tri-diagonal matrix. For the luminosity equation, 


C = -PS^ 
da 


C, + PiS j'^(^ '^ =0. 




and the temperature evolution equation becomes, 




r _ r '-p _ piold 

i?* + = 0 . 


doj 


At 


Thus, the numerical equations to solve are 
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doj+i 


(A.34) 


(A.35) 


(A.36) 

(A.37) 


In this scheme, the unknowns are (7i, Ti, Ti, •'' ■, En-i-,Tn) and the final equations to solve are 
P 2 N -1 instead of F 2 N since we don’t have Tat+i as unknown. We also use the same boundary 
condition as in Penta-diagonal scheme, £0 = 0 and Cn = e^^A tiP? a bT( = e^'i’47ri?^cjB/(77v) • 
Therefore, 


Fi 


F 2 N -1 
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At 


(A.38) 

(A.39) 
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Newton-Raphson iteration method gives the equations, 


(A.40) 

(A.41) 


where 


F2i-i + + B2^-l5T^+^ + C2i-iSC^^^ = 0, 

F2i + A2^5T^^^ + B2i5C^+^ + C2i5Tr+f = 0 , 
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Special case is needed for the boundary grid points. 


(A.42) 

(A.43) 

(A.44) 
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(A.46) 

(A.47) 
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(A.48) 

(A.49) 

(A.50) 


The tri-diagonal matrix becomes 

(Bi Cl 0 • • • 

A 2 82 C 2 0 

0 A 3 Bs C 3 0 


0 • • • 0 A2N-3 B2N-3 C2N-3 

0 • • • 0 A2N-2 B2N-2 

\ 0 • • • 0 A2Ar-i 
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/ 67^+^ \ 
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F 2 N -3 

C 2 N -2 




F 2 N -2 

B 2 N -1 / 




\ +iV-l / 


(A.51) 


Appendix A.3. Comparison 

Each numerical solution (tri-diagonal, penta-diagonal, and tri-diagoal all grids) gives the similar 
solution if the initial condition is identical for each simulation. In the point of view of numerical 


40 




Log(t/yr) 



Log(t/yr) 


Figure A. 18: (Color online) Curves for each numerical method. SkI4 is used to simulate neutron 
star cooling. All three methods shows the identical results. The left figure shows the normal non 
pairing phase and the right hgure shows the superconducting phase. 


cost, tri-diagonal even (Lr)-odd (T) method is superior to penta-diagonal even (Lr)-odd (T) and 
tri-diagoal all grids method. Figure A. 18 shows neutron star cooling curves with SkI4 model. Three 
different numerical methods show almost identical results. The difference in early stage is caused by 
the difference in the time step At in each simulation. That is, penta-diagonal scheme, for example, 
for some case, tri+i /2 is normal state and can be superfluidic phase because of temperature 
difference in each step. Thus the time step should be adjusted to solve the diffusion equations. For 
normal phase, all three methods give no difficulty in the simulation. However, in superconducting 
phase, the most stable numerical method is tri-diagonal with even (L^) and odd (T) scheme since 
it is free from the intermediate time step for sudden decrease of temperature. 


Appendix B. Spatial zone and time step 

In neutron star cooling simulation, we make grids from the core to outer boundary of crust 
{p = lO^^g/cm^) and connect the temperature Tf, with Tg using uniform luminosity approximation 
and Tg — Tb relation IMIETIEI]. The density of the core is around p ~ 10^^ ~ lO^^g/cm^ and the 
crust has the density in the range of 10^*^ to 10^'^g/cm^. Even though, the size of crust is only ~ 
1km, the nuclear phase changes from heavy nuclei with neutron and electron gas to heavy nuclei 
with electron gas. Since the different equation of state gives different central, core-crust boundary, 
and neutron drip density, it is reasonable to make mesh point, 

iVi = VFilogiof—) , (B.l) 

N2 = VFalogiof—) , (B.2) 

As = VFslogio , (B.3) 

\ Penv J 

where pc is the central density, pcore is the density for core-crust boundary, pdrip is the neutron drip 
denisty, and penv = 10^°g/cm^ for density of boundary of crust and envelope. 
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Figure B.19: (Color online) Different number of zone for core and crust region. Large scale figure 
(left) shows that three curves give the same behavior. Enlarge figure near the critical temperature 
shows that W = 100 is enough for the numerical simulation. 


Several constraints for time step At are used. In general, as time goes, the numerical solution 
is more stabilized so that we can use large time step for the numerical solution. Here we use the 
Tgff to determine the next time step. We choose different tscale-, At^^^ = tgcaie^t^ for different 
conditions of T^ff and 


T n rpn—1 
off -L t 


-'ey/ 

> 0.1, 

rj-iTl — 1 

0.05 < 

rpn rpn — l 

-'ey/ 


^^Tl —1 
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pTl— 1 

Cff 

rpn rpn—i. 

< 0.01 


< 

< 


0 . 1 , 

0.05, 


tscale — 1.02 , 
t scale 1.1 , 

tscale — 1.2 , 

^ scale — 1.5 . 


(B.4) 


Another constraint for the time step comes from total time. In our simulation the next time step 
is always less than one tenth of total time, 

Ar+^ = min(tscaieAr, ^t). (B.5) 

If the superfluidity happens, a neutron star experiences drastic changes in specific heat, thermal 
conductivity, and neutrino emission rate. Thus, when the internal temperature drops below the 
critical temperature for superfluidity, we use adaptive time step method. The should change 
within maximum 5% of For instance, if the numerical solution gives < 0.95 we 

solve the diffusion equations again with the new time step (where index 

i indicates the ith trial time step) until > 0.957^y. In our simulation Ueduce = 0.75 to reduce 

the time step. Once we find the solution, according to the temperature differences between and 
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Figure B.20: (Color online) Left figure shows the large scale cooling curve. Right figure shows the 
cooling curve near the critical temperature for superfluidity. Each curve shows different behavior 
near the critical temperature. 


we use the adaptive tscale for the next time For superfluidity case, we use 

<0.05, ts^ = 1.2, 


0.01 < 




rpn _ 'JM 


£// 


ff 


ff 


a—1 

eff 


=// 

< 0.01 


(B.6) 


ts2 = 1.5. 


In Fig. B.20[ we compare results from three different numerical methods. If the superfluidity 
happens, penta-diagonal method needs a smaller size of time step to make the result similar with 
the ones from both the tri-diagonal methods (even Lr and odd T) and the tri-diagonal methods in 
which Lr and T are defined in all grid points. 
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